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LONG-TERM  GOALS 

Our  goal  is  to  provide  a  complete  simulation  code  that  will  represent  and  predict  the  sediment  transport 
and  bed  features  on  the  continental  shelf  at  user-specified  resolution  by  using  state-of-the-art 
algorithms  for  the  physics  and  numerics  of  the  simulation  code. 

OBJECTIVES 

Our  primary  objective  is  to  simulate  the  ripple  climate  on  the  bed  of  the  inner  shelf  at  depths  on  the 
order  of  20  m  and  over  domains  ranging  from  centimeters  to  kilometers  in  support  of  the  Ripples  DRI 
experiments  and  analyses.  Our  secondary  objective  is  to  use  simulation  to  better  understand  the 
physics  of  ripple  formation  and  sediment  transport  in  this  environment.  A  third  objective  is  to  create  a 
powerful  and  fast  parallel  numerical  code  system  by  the  synthesis  of  a  number  of  our  computational 
tools. 

APPROACH 

We  are  developing  a  multiscale  simulation  code  that  will  represent  and  predict  the  sediment  transport 
and  bed  features  on  the  continental  shelf  via  an  integrated  modeling  framework  that  couples  a 
laboratory-scale  large-eddy  simulation  code  with  our  field-scale  nonhydro  static  coastal  ocean  model 
SUNTANS  (Eringer  et  ah,  2006).  Ultimately,  the  wave  climate  as  predicted  by  SUNTANS  will  drive 
the  formation  of  ripples  in  the  small-scale  code,  which  would  in  turn  alter  the  bottom  roughness  as 
seen  by  the  waves,  and  alter  their  behavior  accordingly.  Using  this  tool  we  will  have  the  ability  to 
predict  the  evolution  of  waves  and  benthic  ripples  on  coastal  shelves  and  in  shallow  embayments  given 
the  distribution  of  sediment  properties  and  the  incoming  wave  field.  The  components  of  the  multiscale 
simulation  code  are  described  below. 

Eor  the  field-scale  solver,  we  are  using  the  parallel  unstructured  grid  code,  known  as  SUNTANS 
(Stanford  Unstructured  Nonhydrostatic  Terrain-following  Adaptive  Navier-Stokes  Simulator),  which 
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solves  the  nonhydrostatic  Navier-Stokes  equations  under  the  Boussinesq  approximation  with  the  large- 
eddy  formulation  for  the  resolved  motions.  The  code  has  been  developed  by  Fringer  et  al.  (2006)  in 
order  to  determine  the  behavior  of  the  nonhydrostatic  internal  wave  spectrum  in  Monterey  Bay  (Jachec 
et  al.,  2006).  SUNTANS  is  a  high-resolution  simulation  tool  for  coastal,  estuarine,  and  river 
simulations,  and  incorporates  a  host  of  processes,  from  the  free-surface  to  the  bed,  that  play  an 
important  role  in  driving  smaller-scale  motions,  including  the  nonhydrostatic  pressure,  wind  stress, 
surface  and  internal  tides,  and  parameterizations  of  small-scale  internal  as  well  as  surface  breaking 
physics.  Most  importantly,  because  it  is  unstructured,  and  adaptive  in  the  long-term,  we  have  the 
ability  to  employ  higher  resolution  around  the  site  of  interest  for  the  Ripples  DRI  field  experiments 
where  we  embed  the  smaller-scale  simulation  codes  that  we  describe  below. 

The  SUNTANS  code  is  being  used  to  determine  the  field-scale  motions  that  drive  the  smaller-scale 
flow,  which  is  computed  with  a  nonhydrostatic  parallel  large-eddy  simulation  Navier-Stokes  solver. 
This  solver  employs  the  code  of  Cui  and  Street  (2001;  2004)  (henceforth  referred  to  as  PCUI)  which  is 
an  adaptation  of  the  large-eddy  simulation  solver  of  Zang  et  al.  (1994)  to  simulate  laboratory-scale 
rotating  stratified  flows  using  MPI,  the  message  passing  interface,  on  parallel  computers.  The 
formulation  of  Zang  et  al.  employs  a  non-orthogonal  curvilinear,  boundary-following  grid  as  well  as  a 
dynamic-mixed  subfilter  scale  model  (Zang  et  al.,  1993)  to  compute  the  subfilter  scale  terms  that  arise 
from  volume  filtering.  The  equations  are  discretized  in  time  with  a  semi-implicit,  second-order 
accurate  method  and  in  space  with  second-order  differencing  and  advanced  with  a  fractional 
step/projection  method.  The  method  has  been  applied  to  a  host  of  simulations,  including  turbulent 
stratified  flow  over  a  wavy  bed  (Calhoun  et  al.,  1999),  interfacial  breaking  waves  (Fringer  &  Street, 
2003),  and  the  study  of  a  round  jet  in  crossflow  (Yuan  et  al.,  1999)  using  a  total  of  20  overlapping 
grids. 

For  sediment  transport,  we  are  employing  algorithms  from  the  code  of  Zedler  and  Street  (2001;  2002; 
2006),  who  extended  the  formulation  of  Zang  et  al.  (1994)  and  added  a  module  which  solves  the 
advection  diffusion  equation  (with  settling  term)  for  the  sediment.  Under  this  work  Zedler  and  Street 
also  created  a  bed-load  module  and  a  bed  erosion/deposition  module  (Zedler,  personal  communication, 
2005)  to  go  with  the  module  for  suspended  sediment  transport.  Zedler  provided  us  with  the  equations 
in  the  useful  terrain-following  coordinates  of  her  original  code  so  that  the  resolution  of  shear  stress 
direction,  for  example,  is  clear.  Motion  of  the  bed  in  SUNTANS  will  be  computed  using  a  boundary- 
fitted  curvilinear  grid  that  moves  with  the  bed,  following  the  formulation  of  Hodges  and  Street  (1999). 

WORK  COMPLETED 

We  have  implemented  the  suspended  sediment  transport  model  into  the  PCUI  code,  including  a  new 
bottom  boundary  condition  for  sediment  transport  in  large-eddy  simulation.  We  have  also 
implemented  a  boundary-fitted,  moving  curvilinear  grid  that  follows  arbitrary  geometries  and  enables 
simulation  of  ripple  formation  and  evolution  while  guaranteeing  conservation  of  suspended  sediment 
mass. 

RESULTS 

Our  work  to  date  has  focused  on  incorporating  the  suspended-sediment  modules  of  Zedler  and  Street 
(2001,  2006)  into  the  parallel  large-eddy  simulation  code  PCUI  (Cui  and  Street,  2001)  in  order  to 
perform  high-resolution  simulations  of  sediment-ripple  dynamics  on  parallel  computers.  We  have 
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employed  several  modifications  to  the  implementation  of  Zedler  and  Street,  most  notably  the 
implementation  of  the  bottom  boundary  condition  for  sediment,  which  does  not  require  knowledge  of 
the  near-bed  turbulence  parameterization.  In  the  formulation  of  Zedler  and  Street,  the  bottom 
boundary  condition  for  the  suspended  sediment  is  given  by 


Vj  dC  _ 
Pr  dz 


(1) 


where  Vj.  is  the  eddy-viscosity,  Prj  is  the  turbulent  Prandtl  number,  C  is  the  grid-filtered  sediment 
concentration  (in  a  large-eddy  simulation  framework),  z  represents  the  vertical  direction,  and  Pk  is  the 
pick-up  function.  Although  the  same  formula  was  employed  by  van  Rijn  (1986)  with  a  parabolic  eddy- 
viscosity  model,  the  drawback  to  using  this  boundary  condition  is  that  it  can  be  highly  restrictive  when 
the  magnitude  of  the  eddy-viscosity  is  small,  since,  given  a  value  of  the  pick-up  function  Pt,  the  near¬ 
bed  gradient  of  the  suspended  sediment  concentration  is  given  by 


dC 

dz 


(2) 


Use  of  this  formula  to  specify  the  bottom  boundary  condition  for  the  suspended  sediment 
concentration  is  prohibitive  in  a  large-eddy  simulation  because  the  near-bed  eddy-viscosity  is  small, 
particularly  when  high-resolution  is  employed. 
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Figure  1:  Two-dimensional  schematic  of  a  near-bed  control  volume  of  with  Ax  and  height  Az, 
showing  the  filtered  fluxes  that  influence  the  rate  of  change  of  the  suspended  sediment 

concentration  within  the  control  volume. 


Our  breakthrough  was  achieved  by  considering  a  finite  volume  form  of  the  filtered  transport  equation 
for  the  suspended  sediment.  In  what  follows  we  derive  it  in  two  dimensions  for  illustrative  purposes 
although  we  employ  the  three-dimensional  case.  The  governing  equation  is  obtained  by  filtering  the 
transport  equation  for  the  suspended  sediment  and  discretizing  it  over  a  two-dimensional  finite  volume 
of  width  Ax  and  height  Az,  with  the  east-and  west-face  fluxes  denoted  by  e  and  w,  respectively,  and  the 
top  and  bottom-face  fluxes  denoted  by  t  and  b,  respectively  (see  Figure  1),  to  obtain 
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dC  ^  1 
dt  Ax 


uC 


-uC 


+  - 


Az 


(w-wJC  -(w-wJC 


=  0 


(3) 


where  is  the  settling  veloeity.  The  standard  proeedure  in  a  large-eddy  simulation  is  to  employ  a 

model  for  the  unelosed  terms  uC  and  wC ,  whieh  arise  from  filtering  the  nonlinear  transport  terms.  In 
our  work  we  employ  the  dynamie-mixed  model  of  Zang  and  Street  (1993),  whieh,  for  the  vertieal  flux 
of  sediment,  gives 


wC  =  wC  +  Xi 


(4) 


where  the  vertieal  subfilter-seale  flux  is  given  by  the  sum  of  a  gradient-diffusion  and  the  modified 
Leonard  term  =wC  -  wC  ,  sueh  that 


Z3 


=  wC  -wC  =  -kj 


dz 


+  P". 


(5) 


where  kj.  =  Vj.  /Pr^-is  the  sealar  diffusivity.  Employing  this  model  for  the  unclosed  terms  in  equation 
(3)  requires  that  all  terms  in  equation  (5)  be  specified  at  the  bed  (the  b  term  in  equation  3),  thereby 
requiring  a  gradient-diffusion  type  model  of  the  form  given  in  equation  (2).  In  order  to  avoid  the 
aforementioned  limitation  associated  with  a  small  eddy-viscosity  at  the  bed,  we  specify  a  boundary 
condition  at  the  bed  in  equation  (3)  for  the  term  at  the  bottom  face  before  imposing  a  model  to  handle 
the  unclosed  terms.  Specifically,  in  our  formulation  we  impose  a  flux  on  the  unclosed  term  at  the  bed 
of  the  form 


(w-wJC 


=  E-D 


(6) 


where  the  erosion  is  modeled  with  the  pickup  function  E  =  P,^,  and  the  deposition  is  given  by 
D  =  ,  where  is  the  value  of  the  suspended  sediment  concentration  at  the  bed  that  is 

obtained  via  extrapolation  from  interior  points. 

We  have  employed  this  boundary  condition  to  simulate  the  dynamics  of  suspended  sediment  in  a 
channel  flow  of  length  L=2.0  m,  width  W=1.0  m,  and  height  H=0.6  m,  and  mean  streamwise  velocity 
of  U=1  m  s'  ,  to  yield  a  Reynolds  number  of  600  000  (based  on  Re=UH/v,  where  v  =10'  m  s'  is  the 
kinematic  viscosity  of  water).  The  flow  is  started  from  rest  and  allowed  to  come  to  a  statistically 
steady  state,  whereupon  the  sediment  is  initialized  with  a  uniform  distribution.  As  depicted  in  Figures 
2(a)-(c),  despite  the  existence  of  turbulent  fluxes  at  the  start  of  the  simulation,  and  because  there  is  no 
sediment  supply  at  the  top  surface,  gravitational  settling  dominates  the  behavior  of  the  suspended 
sediment  profiles.  As  a  result,  there  is  a  strong  downward  sediment  flux  and  the  total  concentration 
decreases  rapidly  during  this  initial  settling  period.  Without  turbulent  diffusion,  this  process  would 
persist  until  all  sediment  particles  settle  out.  However,  once  turbulent  transport  becomes  strong 
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enough,  particularly  when  the  near-bed  concentrations  become  high  due  to  sediment  pickup,  sediment 
is  transported  from  the  channel  bed  into  the  upper  part  of  the  water  column.  As  depicted  in  Figures 
2(c)-(d),  turbulent  diffusion  becomes  significant  and  this  results  in  a  net  upward  flux  of  sediment. 
During  this  transitional  period,  due  to  turbulent  diffusion,  the  concentration  field  gradually  increases 
with  time  due  to  erosion  until  the  upward  turbulent  flux  balances  the  downward  flux  due  to  settling.  At 
this  point,  the  suspended  sediment  concentration  reaches  a  statistically  steady  state  and  profiles  at  each 
cross  section  remain  statistically  invariant  (see  Figures  2(e)-(f)). 


Figure  2:  The  temporal  evolution  of  the  laterally-averaged  suspended  sediment  concentration 
profiles,  normalized  by  the  measured  reference  concentration  Ca  at  0. 02  m  above  the  bed.  Each  line 
represents  a  different  streamwise  location  along  the  channel. 


The  strong  spatial  variability  of  the  suspended  sediment  concentration  is  shown  in  Figure  3,  which 
depicts  the  concentration  contours  in  a  vertical  plane  along  the  centerline  at  t  =  1002  s.  In  order  to 
further  investigate  how  sediment  suspension  is  linked  to  the  ambient  flow  field,  in  Figure  4  we  plot 
instantaneous  concentration  contours  superimposed  over  the  velocity  vectors  in  a  cross-sectional  plane 
at  X  =  0.625  m.  The  figure  shows  transport  of  the  suspended  sediment  associated  with  the  resolved 
large  eddies.  Near  the  bottom,  patches  of  high  sediment  concentration  correspond  to  high  eddy 
intensity  or  upward  flux.  These  eddies  form  vortex  cores  and  can  extend  downstream,  thereby  greatly 
enhancing  sediment  transport  in  the  streamwise  direction.  Because  the  mean  vertical  sediment  flux  due 
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to  the  eddy-diffusion  is  mueh  less  than  that  induced  by  the  resolved  field,  the  transport  of  suspended 
sediment  in  our  results  is  mainly  due  to  advection  by  the  resolved  eddies.  Once  sediment  is  mobilized 
by  the  excess  shear  stress,  which  we  model  with  the  pick-up  function,  it  is  further  suspended  into  the 
water  column  only  if  the  ambient  turbulent  flux  is  strong  enough  to  overcome  the  gravitational  settling. 
Otherwise,  suspended  sediment  will  deposit  at  the  bottom  and  its  motion  will  be  dominated  by 
granular  forces.  This  eddy-transport  phenomenon  is  illustrated  further  by  a  zoom-in  plot  of  Figure  4  in 
Figure  5,  which  depicts  high  sediment  concentration  with  weak  spatial  variability  near  the  bed  and 
suspension  of  sediment  due  to  the  strong,  coherent  eddies  and  upward  velocity.  These  near-bottom 
vortex  structures  and  the  associated  sediment  transport  are  also  found  by  Zedler  and  Street  (2001, 
2006),  who  simulate  high  Reynolds  number  channel  flow  using  LES  with  a  different  wall  model  and 
boundary  conditions. 
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Figure  3:  A  snapshot  of  the  concentration  contours  in  a 
vertical  plane  along  the  channel  centerline  at  t  =  1002  s. 
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Figure  4:  Velocity  vectors  superimposed  over  the  concentration 
contours  in  a  cross-sectional  plane  at  x  =  0.625  m  and  t  =  1002  s. 
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Figure  5:  A  zoom-in  view  of  Figure  4  in  the  region  where  high  concentration 
is  associated  with  the  strong  near-bottom  eddy  and  upward  flow. 
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